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ESTIMATION OF THE GRADIENT DENSITY IN ARBITRARY 
FINITE DIMENSIONS USING THE STATIONARY PHASE 
APPROXIMATION 

KARTHIK S. GURUMOORTHY**, JOHN CORRING* 1 ", AND ANAND RANGARAJAN" 

Abstract. We prove a novel result wherein the density function of the gradients of a smooth 
function S (obtained via a random variable transformation of a uniformly distributed random vari- 
able) defined on a compact domain Q C WL d is accurately approximated by the normalized power 

^ 1 spectrum of <f> = exp [rfj as the f ree parameter r — > 0. The result is shown using the well known 

_ stationary phase approximation and standard integration techniques and requires proper ordering of 

limits. Experimental results provide anecdotal visual evidence corroborating the result. 
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1. Introduction. 

| The centerpiece of our current work is to provide a useful application of the 

stationary phase method to compute the joint density function corresponding to the 
gradients of the smooth function S (joint density function of VS). Here, the joint 
density function is obtained via a random variable transformation of a uniformly 
^ | distributed random variable X over the connected compact domain il C R d using the 

gradients VS* as the transformation function. In other words, if we define a random 
variable Y = VS(X) where the random variable X has a uniform distribution on f2, 
the density function of Y represents the joint density function of the gradients. 
QO | In this paper we show that the density of the gradients of S can be computed 

C^) ■ indirectly, without explicitly determining VS*, by expressing S as the phase of a wave 

function </>, specifically </>(x) = exp ^ ^ fo r small values of t, and then consider 
the normalized power spectrum — magnitude squared of the Fourier transform — of 4> 
[3J. Using the stationary phase approximation — a well known technique in asymptotic 
\ analysis |10| — we show that in the limiting case as r — > 0, the power spectrum of 4> 

converges to the density of Y and hence can serve as its density estimator at small, 
non-zero values of r. In other words, if P(u) denotes the density of Y and if -P-r(u) 
corresponds to the power spectrum of <\> at a given value of r, Theorem 13 . 2 1 constitutes 
the following relation namely, 



lim / P T (u)du= / P(u)du 



(1.1) 



where A/^(uo) is a small neighborhood around uo. We call our approach as the 
wave function method for computing the probability density function and henceforth 
will refer to it by the same name. We would like to emphasize that our work is 
fundamentally different from estimating the gradients of a density function [S] and 
should not be semantically confused with it. 
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1.1. Motivation from quantum mechanics. Our new mathematical rela- 
tionship is motivated by the classical-quantum relation, wherein classical physics is 
expressed as a limiting case of quantum mechanics [7J H] . When S is treated as the 
Hamilton- Jacobi scalar field, the gradients of S correspond to the classical momentum 
of a particle [6]. In the parlance of quantum mechanics, the squared magnitude of 
the wave function expressed either in its position or momentum basis corresponds to 
its position or momentum density respectively. Since these representations (either 
in the position or momentum basis) are simply (suitably scaled) Fourier transforms 
of each other, the squared magnitude of the Fourier transform of the wave function 
expressed in its position basis is its quantum momentum density. However, the time 
independent Schrodinger wave function 0(x) (expressed in its position basis) can be 

approximated by exp ( as t -> ^. Here r (treated as a free parameter in 

our work) represents Planck's constant. Hence the squared magnitude of the Fourier 

transform exp ^ ^ corresponds to the quantum momentum density of S. The 

principal results proved in the article (Theorem I3.2[) state that the classical momen- 
tum density (denoted by P) can be expressed as a limiting case (as r — > 0) of its 
corresponding quantum momentum density (denoted by P T ), in agreement with the 
correspondence principle. 

2. Existence of joint densities of smooth function gradients. We begin 
with a compact measurable subset fi of W 1 on which we consider a smooth function 
S : f2 — > R with derivatives of arbitrary order. We use /i to denote Lebesgue measure 
on M. d . We assume that S is well-behaved on the boundary of fl as elucidated in 
Appendix [B] Let H x denote the Hessian of S at a location x <E il and let det (H x ) 
denote its determinant. Let cr x correspond to the signature of the Hessian of S at x 
defined as the difference between the number of positive and negative eigen-values of 
"H x - In order to exactly determine the set of locations where the joint density of the 
gradients of S exists, consider the following three sets: 



A u = {x : V5(x) = u} (2.1) 
B = {x : det(Hx) = 0} (2.2) 
C = {WS(x) :xe BUdn}. (2.3) 

Let iV(u) = |.4 U |. We employ a number of useful lemma, stated here and proved 
in Appendix [A"l 

Lemma 2.1. [Finiteness Lemma] A u is finite for every u ^ C. 

Lemma 2.2. [Neighborhood Lemma] For every Uo ^ C, there exist a closed neigh- 
borhood M v (uq) around Uq such that Af v (uo) C\C is empty. Furthermore, i/|_4 u0 | > 0, 
we can find a closed neighborhood A/^(x) around each x £ _4 u0 such that 

1. VS (XM) = M^oj 

2. det (H y )^0,VyeJV,(x) 

3. The inverse function V '5 X 1 '(u) : A/"^(uo) — > J\f n {x) is well-defined. 

4. For y, z G J\f v (x),a y = er z . 
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Lemma 2.3. [Density Lemma] Given S G C°°(fi), X ~ UNI(Q), the probability 
density of Y = WS(X) on M. d — C is given by 

1 N(u) 1 

p(u > ^ MH) g P^mtO (2 - 4) 

w/iere x fc e _4 u ,Vfc G {1,2, • • • , iV(u)}. 



From Lemma 12.31 it is clear that the existence of the density function P at a 
location u G M. d necessitates a non-vanishing Hessian matrix (det('H) ^ 0) Vx <E .4 U . 
Since we are interested in the case where the density exists almost everywhere on 
M. d , we impose the constraint that the set B in Equation 12.21 comprising of all points 
where the Hessian vanishes, has a Lebesgue measure zero. It follows that /u(C) = 0. 

2.1. Characteristic function formulation for obtaining densities. The 

characteristic function ^y(a>) for the random variable Y = VS(X) is defined as the 
expected value of exp (iuj.Y), namely 

i> Y («) = E [exp {ioJ.Y)] = -4yT [ exp (ico ■ VS(x)) dx. (2.5) 

M^') Jn 

Here denotes the density of the uniformly distributed random variable X on fi. 

The inverse Fourier transform of characteristic functions also serves as the density 
functions of random variables [I]. In other words, the density function P(u) of the 
random variable Y can be obtained via 



P (u) = t^w y Vv(w)exp(-£w • u)dw 

exp {iu ■ [V5(x) - u]} dxdw. (2.6) 



(2tt) 
1 

~ (27r)«V(n) y 7n 

In this work we showcase a direct relation between the characteristic function formu- 
lation for computing densities and the power spectrum of the wave function </>(x) = 

exp ^ tS | x ' ) ^ at small values of r, as elucidated in Section 13721 Nonetheless, the sta- 
tionary phase method has little bearing with the characteristic function approach for 
estimating the probability densities and should be treated as an independent recipe 
for computing the density function. Though we build an informal bridge between the 
wave function and the characteristic function approaches in Section 13. 2\ based on the 
reasons described therein, we strongly believe that the mechanism of stationary phase 
approximation is essential to formally prove our main Theorem 13.21 Hence our wave 
function method should not be merely treated as a reformulation of the characteristic 
function approach. 

3. Equivalence of the Density of Gradients and the Power Spectrum. 

Define the function F T : M. d — >• C as 

Mu) = J^hsW* In eXP {~r [5(X) " U ' X] } ^ (3 ' 1} 

for t > 0. F T is very similar to the Fourier transform of the function exp{ tlS | x ) }, 
The normalizing factor in F T comes from the following Lemma whose proof is given 
in Appendix \X\ 
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Lemma 3.1. [Integral Lemma ] F T e L 2 (R d ) and \\F T \\ 2 = 1. 
The power spectrum is defined as [2] 



P T (u) = P r (u)F r (u). (3.2) 

Note that P T > 0. From Lemma (|3.1j) we see that / P r (u)du = 1. Hence, treating 
P T {u) as a density function, we have the following theorem. 
Theorem 3.2. For u ^ C, 

lim ,.A ^ lim /" P T (u)du = P(u ) (3.3) 
«->o ^ (7V Q (u )) r->o y 

A/a(ll ) 

where A/" a (uo) is a small ball around Uo o/ radius a. 

Before embarking on the proof, we would like to emphasize that the ordering of the 
limits and the integral as given in the theorem statement is crucial and cannot be 
arbitrarily interchanged. To press this point home, we show below that after solving 
for P T , the lim r _yo Pt does not exist. Hence, the order of the integral followed by the 
limit r — > cannot be interchanged. Furthermore, when we swap the limits between 
a and r, we get 

lim lim — ttt—, — rr f P T (u)du = lim P T (u Q ) (3.4) 

r _ >0 a->0/i(A/a(u )) J r^O 
JV Q (u ) 

which does not exist. Hence, the theorem statement is valid only for the specified 
sequence of limits and the integral. 

3.1. Brief exposition of the result. To understand the result in simpler terms, 
let us reconsider the definition of the scaled Fourier transform given in Equation [3TTJ 
The first exponential exp ( tS ^ \ is a varying complex "sinusoid", whereas the second 

exponential exp ( x ) is a fixed complex sinusoid at frequency —. When we multiply 
these two complex exponentials, at low values of r, the two sinusoids are usually not 
"in sync" and tend to cancel each other out. However, around the locations where 
VS'(x) = u, the two sinusoids are in perfect sync (as the combined exponent is sta- 
tionary) with the approximate duration of this resonance depending on "H x . The value 
of the integral in (|3.ip can be approximated via the stationary phase approximation 
[TO] as 

N(u) ( . 



F ^ * h exp i - [5(Xfc) " u ' Xfc] + mxk ^ / 7mm\ - (3 - 5) 

The approximation is increasingly tight as r — ¥ 0. The squared Fourier transform (P T ) 
gives us the required result i det(n — JT exce Pt f° r the cross phase factors 

5(xfc) — 5(x;) — u • (xfe — x;) obtained as a byproduct of two or more remote locations 
Xfe and xj indexing into the same frequency bin u, i.e, x^ ^ x;, but VS'(xfc) = 
VS*(x;) = u. Integrating the squared Fourier transform over a small neighborhood 
A/" Q (u) around u removes these cross phase factors and we obtain the desired result. 
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3.2. Relation between the characteristic function and the power spec- 
trum formulation of the density of gradients. We showcase the relationship 
between the commonly used characteristic function formulation of the density and 
our formulation arising from the power spectrum. We request the reader to bear in 
mind the following precautionary note. What we show below cannot be treated as a 
formal proof of Theorem 13.21 The main motivation of this section is to provide an 
intuitive reason behind our theorem, where we try to directly manipulate the power 
spectrum of </>(x) into the characteristic function formulation given in Equation 12.61 
circumventing the need for the closed-form expression of the density function P(u) 
given in Equation 12.41 Though we do not obtain a formal proof through this ap- 
proach, our attempt provides a mathematically intuitive approach to establishing the 
equivalence between the power spectrum and the characteristic function formulation 
and thereby to the density function P{u). The formal proof is given in the subsequent 
section by using the method of stationary phase approximation. 

For simplicity we choose to consider a region fi that is the product of closed 

d 

intervals, fl = J\ [o-i,bi]. Based on the expression for the scaled Fourier transform 
F T (u) in Equation 13.11 the power spectrum P T (u) is given by 



P ^ (U) = (2^km L L ^ { ^ [S{X) " S(y) - (U - { *- y))1 



Define the following change of variables namely, 



x-y x + y 
uj = , v = — - — . 



(3.7) 



Then, 



UJT UJT 

x = v + — ,y = v 



2 2 
and the integral limits for uj and v are given by 



(3.8) 



i=l 
i! 

K,=n 



a, — hi hi - a, 



i=l 



(3.9) 
(3.10) 



where uii is the i th component of uj. Note that the Jacobian of this transformation is 
r d . Now we may write the integral P T (u) in terms of these new variables as 



P T (u) 



(27r) d /i(ft) 



(3.11) 



where, 

£(w,u) 



exp 



S 



("+t) 



S u 



cxp {— iu ■ u>} dh>. (3.12) 
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The mean value theorem applied to 5 (1/ + ^ ) ~ S (is — ^) yields that 

f(«,u)= fexp{iu- [VS(z(u,u>))-u]}du (3.13) 

where, 

z( VjW ) = c(i/+^) + (l_ c )( I/ _Z^) (3.14) 

with c € [0, 1]. When u> is fixed and r — > 0, z(i/, <*;)—>• 1/ and so for small values of r 
we get 

£(w,u)w y exp{iw[V5(i/)-u]}di/. (3.15) 

Again we would like to drive the following point home. We do not claim that we have 
formally proven the above approximation. On the contrary we believe that it might 
be an onerous task to do so as the mean value theorem point z in Equation 13.131 is 
unknown and the integration limits for v directly depend on r. The approximation is 
stated with the sole purpose of providing an intuitive reason for our Theorem 13.21 

But note that the integral range for u> depends on r and so when u = O (i), 
ujt -ft as t —> and hence the above approximation for £(oj,u) in Equation 13.151 
might seem to break down. To evade this seemingly ominous problem, we manipulate 
the domain of integration for uj as follows. Fix an e G (0, 1) and let 



W = W°° U W T = yW \ Y[[-Mi, Mi] ) U JJ[-M is Mi] (3.16) 

where, 



i=i 



M i = (b i -a i )T £ - 1 . (3.17) 

Observe that in W T , uj is deliberately made to be 0(r e_1 ) and hence u>t — >• as r — > 0. 
Hence the approximation for ^(uj,u) in Equation 13.151 might hold for this integral 
range of lj. For consideration of u> € W°°, recall that the Theorem 13.21 requires the 
power spectrum P T (u) to be integrated over a small neighborhood Af a (uo) around 
uo- By using the true expression for u) from Equation 13. 121 and performing the 
integral for u prior to u) and u 1 we get 

£(u>, u)dujdu = 



JV„(uo) W" 



J j exp|iw- V5 (c(u + ^) + (1 — c)(u — -^-)) } J ex P { — ^ ' u } dudvdu; 

(3.18) 

Since both the M, in Eq uation 13 . 1 71 and the lower and the upper limits for u>i, ± bi ~ ai 
respectively, approach 00 as t -> 0, the Riemann-Lebesgue lemma guarantees that 
Mu) 6 W°°, the integral 

J cxp {— iu) ■ u} du (3.19) 

A^c(uo) 
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approaches zero as r — > 0. Hence for small values of r, we can expect the integral over 
W T to dominate over the other. This leads to the following approximation namely, 

A/- Q (u ) A(u )W" 

as t approaches zero. Combining the above approximation with the approximation 
for £(u>,u) given in Equation 13.151 and noting that the integral domain for u> and 
v approaches M. d and fl respectively as r — > 0, the integral of the power spectrum 
P T {u) over the neighborhood A/" Q (uo) at small values of r in Equation 13.111 can be 
approximated by 

j P T {u)du^ (2?r) ^ (n) J JJexp{iu>- [VS(is)-u]}dvdudu. (3.21) 



This form exactly coincides with the expression given in Equation l2.6l obtained through 
the characteristic function formulation. 

We believe that the approximations given in Equations 13.151 and 13.201 cannot 
be proven easily as they involve limits of integration which directly depend on r. 
Furthermore, the mean value theorem point z in Equation [3T3] is arbitrary and cannot 
be determined beforehand for a given value of r. The difficulties faced here also 
emphasizes the need for the method of stationary phase approximation to formally 
prove Theorem 13. 2[ as established in the subsequent section. 



3.3. Formal Proof. We wish to compute the integral 

1 



rf /!„,oM/2 / ex P U(5(x)-u-x) 



v. 



at small values of r and exhibit the connection between the power spectrum P r (u) 
and the density function P(u). To this end define ^(x; u) = 5(x) — u ■ x. 



case (i): We first consider the case where N(u) = 0, i.e. u ^ ViS(fi). In other words 
there are no first kind stationary points [10J for this value of u. The proof that this 
case yields the anticipated contribution of zero follows clearly from a straightforward 
technique commonly used in stationary phase expansions when the contribution from 
level sets along the boundary is immaterial. We assume that the function S is suf- 
ficiently well-behaved on the boundary such that the total contribution due to the 
second kind stationary points |10| approaches zero as r — > 0. To keep up with the 
crux of our work and also to improve clarity, we reserve the discussion concerning the 
well-behavior of S and the second kind stationary points in Appendix [B] Contribu- 
tions from the third kind stationary points can also be ignored as they are negligible 
compared to the first kind and approach zero as r tends to zero |1Q[ |S] . 

Lemma 3.3. Fix u i C. If A u „ = then F T (u a ) = 0{y/r) as t -» 0. 

Proof. Recall that there are no first kind stationary points, i.e, V^(x; uo) 7^ 0, Vx. 
Hence the vector field 
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is well-defined. By considering an expression similar to Equation IA.11I — where we 
replace j?fc,z( u ) by ^(x; u ) — and applying the divergence theorem n > d/2 times, the 
Fourier transform in Equation 13.221 can be rewritten as 

Fr(Uo) = (2^W^(OW2 3n(x)cxp|^(x)|dx (3.24) 

Q 

I7a E( iT ) J '/ [v J (x(y))-n]exp|l*(x(y))|dy, (3.25) 



(2 7 rr)^/^(0) 

J — 1 an 



where 



ffj(x) = V • Vj(x) and (3.26) 
v i+i( x ) = |jvfp** (x) - (3 - 27) 

Here n is the unit outward normal to the positively oriented boundary dtt parame- 
terized by y. As n > d/2, the term in the right side of Equation 13.241 is o(^/t) and 
hence can be overlooked. To evaluate the remaining integrals within the summation 
in Equation I3.25[ we should take note that the second- kind stationary points of ^(x) 
on £1 correspond to the first kind stationary points of ^(x(y)) on the boundary dfl 
whose contribution on a d — 1 dimensional space (the dimension of the boundary) is 
0(t^~) as shown in case (ii). After multiplying and dividing by the corresponding 
T "' (j > 1) an d T d / 2 factors respectively, it is easy to see that the total contribution 
of the n integrals in Equation 13. 251 is 0(yfr). Here, we have safely ignored the contri- 
butions of the third-kind stationary points as they are inconsequential compared to 
the other two kinds [10]. Additionally, |F T (uo)| < y/rj(uo) where 7(110) is a bounded 
continuous function of u. For a detailed exposition of the proof, we encourage the 
reader to refer to Chapter 9 in |10| . 
□ 

We then get P T (uo) = 0{t) which can also be bounded by a function of u at small 
values of r. Since V5(f2) is compact and Uo ^ VS'(fi), we can choose a neighborhood 
Af v (uo) around u such that for u € Af v (u ), no first kind stationary points exist and 
hence 

Urn J P T (u)du = 0. (3.28) 

A/^(uo) 

Since N(u) = 0,Vu e A/" l; (u ), the true density P(u) of Y=VS(X) given in Equa- 
tion is also zero for u £ TV^Uo). 

case (ii): Let for uo ^ C, N(uo) > 0. In this case the number of first kind stationary 
points in the interior of Q is non-zero and finite as a consequence of Lemma 12.11 We 
can then rewrite 

1 JV(u) r ( i 1 

F r (u ) = G+ (27rr)rf/2 ^ )1/2 g J exp|-(S(x)-u -x)|dx (3.29) 

fc_1 A^(x fc ) 
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where, 

G = J exp | % - {S{x) - u • x) | dx. (3.30) 

K 

N(u) , ^(u) 

Here the domain K = £1\ 1J Af v (xk). The closed regions < Af n (xfc ) > are obtained 

fc=l J 

from Lemma \2 .21 

By case (i), G = ei(uo,r) = 0{y/r) since K contains no stationary points by 
construction. Also note that the contribution from the overlapping boundaries be- 
tween K and J\f n (xk) cancel each other out as they are orientated in opposite ways. 
Hence we can exclude the contributions from the boundaries of A/^Xfc) — those ob- 
tained from the second and third kind stationary points — while evaluating F t (uq) in 
Equation 13.291 Taking into account the contribution from the first kind stationary 
point at Xfc, we get [10] 

f (i 1 a(?7n-) d / 2 (i 7rl 

_j_ cxp ( r [ * (x) - u ■ x] j dx = vmm\ exp u* (Xfc; Uo) + l(Txt 4 / 

+ e 2 (u ,r) (3.31) 



where £2(uo,t) = O yr 2 J < t 2 72 (uo) f° r a continuous bounded function 72(11). 
Here a G {1/2,1} contingent on whether the first kind stationary point is on the 
boundary or in the interior of f2. Since u ^ C, the first kind stationary points do not 
occur on the boundary and hence a = 1 in our case. Here tr Xfc is the signature of the 
Hessian at Xfc. 

Coupling Equations 13.291 13.301 and 13.311 yields 

1 N{u) f i 11 
^ (U0) = £ exp|^( Xfc;Uo ) +igx ^| ^ det( ^ )| +£3 ( Uo , T ) (3.32) 

where 

e 3 (u , T ) = ei (u ,r) + ( J^ )1/2 = (3.33) 
Based on the definition of the power spectrum P T in Equation 13.21 we get 
1 2V(u) 1 

^(U ) = —TTyT V 1 , (3.34) 
M(^) fr[ |dctCH x J| 

| 1 j ^ ) ^ ) exp{^[vl/(x fc ;uo)-vl/(x ; ;uo)]+z(a Xfc -a Xi )f} 
h h ~ V|dct(H x J|V|det(^ Xi )| 

+ e 4 (u ,r) (3.36) 

where £4(110 includes both the magnitude square of £3(110,1") and the cross terms 
involving the main (first) term in Equation 13.321 and 63(110, r). Notice that the main 
term in Equation 13.321 can be bounded independent of r as 



% IT 

cxp<{ -ty(x k ;u Q ) + ia Xk - 



= l,Vr (3.37) 



10 JOHN CORRING AND KARTHIK S. GURUMOORTHY AND ANAND RANGARAJAN 



and det(H Xk ) ^ 0,Vfc. Since e 3 (u ,r) = 0(y/r), we get e 4 (uo,r) = 0(y/r). Fur- 
thermore, as e 4 (uo,r) can be also be uniformly bounded by a function of u for small 
values of t, we get 

lim / e 4 (u ,T)cZu = 0. (3.38) 



Af„(u ) 

Observe that the Equation l3.34l matches the anticipated expression for the density 
function P(u ) given in Equation l2.41 The cross phase factors in Equation 13.351 arises 
due to multiple remote locations and x/ indexing into u. Since lim r _>o cos (i) is 
not defined, lim T _>.o P t (uq) does not exist and hence the cross cosine terms in Equa- 
tion [3]35]io not vanish when we take the limit. However, the following lemma, which 
invokes the inverse function VSl _1 ^(u) : J\f v (u ) — > A/^(x) defined in Lemma l2~2l 
where x is written as a function of u, provides a simple way to nullify the cross phase 
factors. Note that since each VSi is a bijection, N(u) doesn't vary over Af n (uo). 
Pursuant to Lemma l2.2[ the Hessian signatures cr Xfc ( u ) and cr X; ( u ) also remain constant 
over Af ri (u a ). 

Lemma 3.4. [Cross Factor Nullifier Lemma] The integral of the cross term in 
Eauation \3. 35\ over the closed region A/^(uo) approaches zero as r — > 0, i.e, Vfc 7^ I 

"° J |det(Wx,(»))l 1 ' ! '|det(W Il( „ ) )|V» 

W„(u ) 

The proof is given in Appendix [3] Combining Equations 13.391 and 13.381 yields 

r 1 r N[u) 1 r 

J M(») 7 ^ |det(^ Xfc(u) )| 7 

M,(u ) A/"r,(u ) A/^(u ) 

Letting a < rj completes the proof. 

4. Experimental evidence in 2D. Below, we show comparisons between our 
wave function approach and the standard histogramming technique on some trigono- 
metric and exponential functions sampled on a regular grid between [—0.125, 0.125] x 
[—0.125,0.125] at a grid spacing of 2T3. For the sake of convenience, we normalized 
the functions such that its maximum gradient magnitude value (HVS 1 !!) is 1. Us- 
ing the sampled values S, we computed the Fast Fourier transform of exp (j^j a t 
t = 0.00004, took its magnitude square and then normalized it to compute the joint 
gradient density. We also computed the discrete derivative of S at the grid locations 
to obtain the gradients (S Xl -,S X2 ) and then determined its gradient density by his- 
togramming. For better visualization, we marginalized the density along the radial 
and the orientation directions. The plots shown in Figure (|4.1[) provide anecdotal em- 
pirical evidence corroborating our Theorem 13.21 in 2D. Notice the near-perfect match 
between the gradient densities computed via the standard histogramming and our 
wave function method. 
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Figure 4.1. Comparison results, (i) Left: Gradient magnitude density, (ii) Right : Gradient 
orientation density. In each sub-figure, we plot the density function obtained from histogramming 
and our Fourier transform method in the top and bottom part respectively. 



5. Computational complexity. Though not central to the main theme of our 
current work, we would like to add the following note on the computational time 
complexity of our wave function method in comparison to the characteristic function 
formulation. Given the N sampled values S and its gradients VS, the characteristic 
function defined in Equation 12.51 needs to be computed for A" integral values of u>. 
Each value of uj requires summation over the N sampled values of exp {iu)VS(x)}. 
Hence the total time required to determine the characteristic function is 0(N 2 ). The 
joint density function of the gradients is obtained via the inverse Fourier transform 
of the the characteristic function, which is an 0(N log N) operation. The overall 
time complexity is therefore 0(N 2 ). For our wave function method, the Fourier 

transform of exp j^S^x) j at a given value of r can be computed in 0(N log N) and 
the subsequent squaring operation to obtain the power spectrum can be performed 
in O(N). Hence the density function can be determined in 0(N log AT). Our method 
compares favorably with respect to the characteristic function approach in the overall 
time complexity. Also note that our wave function approach recovers the gradient 
density function P T without explicitly computing the gradients of the function S. 
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Appendix A. Proof of Lemmas. 

1. Proof of Finiteness Lemma 

Proof. We prove the result by contradiction. Observe that A u is a subset of the 
compact set ft. If „4 U is not finite, then by Theorem (2.37) in [9], A u has a limit point 
xo € fi. If xo € dil, then u 6 C giving a contradiction. Else, consider a sequence 
{xnj^j, with each x„ 6 A Ul converging to xo. Since V5 , (x„) = u, Vn, from the 
continuity it follows that VS'(xo) = u and hence xo G *4 U - Let p„ = x„ — Xo and 

h « = Then 

Um V5(x w )-V5(xo)-H XoP » =0 

n^-oo ||p„|| 

where the linear operator "H XQ represents derivative of the function VS : M. d — > M. d at 
the location xo, in other words the Hessian of S at xo. As VS(x„) = VS'(xo) = u 
and H Xo is linear, we get 

lim Hx h n = 0. (A.2) 

n— >oo 

Since h n is defined to be an unit vector, it follows that H Xo is rank deficient and 

det ('Hxo) = 0. Hence x <G B and ueC resulting in a contradiction. 

□ 

2. Proof of Neighborhood Lemma 

Proof. Observe that the set B defined in Equation 12.21 is closed because if xo is a 
limit point of B, from the continuity of the determinant function we have det (T-L XQ ) = 
and hence Xo <E B. Being a bounded subset of it is also compact. As dfl is also 
compact and VS* is continuous, C is compact and hence M. d — C is open. Then for 
Uo ^ C, there exists an open neighborhood Af r (uo) for some r > around Uo such 
that Af r (u ) n C = 0. By letting r\ = |, we get the required closed neighborhood 

A^(u ) C A/" r (u ) containing u . 

Since det (W x ) ^ 0,Vx <G -4 u oi the points 1, 2 and 3 of the lemma follow directly 
from the inverse function theorem. As |.4 u ol is finite by Lemma 12. 1[ the closed 
neighborhood A/^(uo) can be chosen independent of x £ A u q so that the points 1 and 
3 are satisfied Vx € A u q. In order to prove point 4, note that the eigen- values of H x 
are all non-zero and vary continuously for x € A/^x). As the eigen- values never cross 
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the zero mark, they retain their sign and so the signature of the Hessian stays fixed. 



3. Proof of Density Lemma 

Proof. Since the random variable X is assumed to have a uniform distribution on 
il, its density at every location x G Q equals . Recall that the random variable 
Y is obtained via a random variable transformation from X, using the function VS. 
The Jacobian of VS at a location x G equals the Hessian H x of the function S at x. 
Barring the set C corresponding to the image (under V5) of the set of points B where 
the Hessian vanishes and the boundary d£l, the density of Y exists on u G M. d — C and 
is given by Equation 12.41 The reader may refer to pQ for a detailed explanation. 
□ 

4- Proof of Integral Lemma 

Proof. Define a function H(x) by 



Let /(x) = H(x) exp 

Fr(u) 



igfx] 



ff(x) = 

Then, 
1 



1 : if x G O; 
: otherwise. 



/(x) exp 



-i(u • x) 



dx. 



(27TT) d /2^(n) 1 /2 

Letting v = ^ and G(v) = F T (u), we get 

(r) d / 2 M (fi) 1/2 G(v) = | /(x)ex P Hv-x)dx. 

As / is i? 1 integrable, by Parseval's theorem [5] we have 

y |/(x)| 2 dx = r d /i(^) J \F t {vt)\ 2 dv = J \F T (xi)\ 2 dxi. 



By noting that 

/ l/(x)| 2 rfx: 

we get the desired result namely, 



exp 



iS(x) 



g£x = yU(^), 



|F T (u)| 2 du= 1. 



(A.3) 



(A.4) 



(A.5) 



(A.6) 



(A.7) 



5. Proof of Cross Factor Nullifier Lemma 

Proof. Let pk,i{u) denote the phase of the exponential in the cross term barring 
the constant signatures, i.e, 



p fc ,;(u) = $(x fc (u);u) - *(x ; (u);u) 

= 5(x fe (u)) - 5(x ; (u)) - u • (xfe(u) - x m (u)). 



(A- 



Its gradient with respect to u equals 



Vp M (u) = J Xfc [VS(x fc (u)) - u] - J Xi [VS(x,(u)) - u] - x fe (u) + x,(u) (A.9) 
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where J Xk is the Jacobian of x(u) at x k whose (i,j) th equals Similarly for J Xk . 
Since V5(x fc (u)) = VS(xj(u)) = u, we get Vp k ,i(u) = x,(uf- x fc (u) ^ 0. This 
means that the phase function of the exponential in the statement of the lemma is 
non- stationary and hence do not contain any first kind stationary points. Let 

Qk z(u) = - 1 - (A. 10) 

V /|dct(H Xt(u) )| v /|dct(H Xl(u) )| 

Since Vpk.i ^ 0, consider the vector filed pfe i ;(u) = y v^j'ffip 1k,l ( u ) an d note that 
exp \ -p k ,i(u) } q k ,i( u ) = iT I v ' Pfe,/( u )] ex P 1 



irV 



p fej ;(u)exp <^ -pfe,z(u 



(A.ll) 



Inserting Equation IA.11I in Equation 13.351 and applying the divergence theorem we 
get 

J exp |^-p fe ,/(u)| q k j(u)du = it J [V • p k j(u)} exp |^-p fe ,/(u)| alu 



7V,(u ) A^(u ) 



(p fcll ■ n) exp <j -p fei/ (u(v)) } dv. (A.12) 



9Af„(u ) 



Here n is the unit outward normal to the positively oriented boundary cW, ; (uo) pa- 
rameterized by v. In the right side of Equation lA.121 notice that all terms inside the 
integral are bounded. The factor r outside the integral ensures that 

lim j exp j^pfc,;(u)j q kil (u)du = 0. (A.13) 

Af„(u ) 



Appendix B. Well-behaved function on the boundary. One of the foremost 
requirement for the Theorem l3.2l to be valid is that the function ^(x; u) = 5(x) — u • x 
has a finite number of second kind stationary points on the boundary for all-most all 
u. Second kind stationary points are the critical points on the boundary T = dfl where 
a level curve 0(x) = c touches T for some c |101[TT1 [5]. Contributions from the second 

kind are generally O ( T_ ^ _ ) : but an infinite number of them could produce a combined 

effect of O (r d / 2 ), tantamount to the first kind stationary point [TO]. If so, we need 
to account the contribution from the boundary which could in effect invalidate our 
theorem. However, the condition for the infinite occurrence of second kind stationary 
points is so restrictive that for all practical purposes they can be ignored. If the given 
function S is well-behaved on the boundary in the sense explained below, these thorny 
issues could be sidestepped. We would like to emphasize that as we will be integrating 
over u to remove the cross-phase factors, it suffices that the aforementioned finiteness 
condition be satisfied for all-most all u instead of Vu. 
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Let the location x <= T be parametrized by the variable y, i.e, x(y). Let Q(x) 
denote the Jacobian matrix of x(y) whose (i,j) th entry is given by 



Q*j(x) = 




(B.l) 



The second kind stationary points occur at locations x where V^(x(y); u) = which 
translates to 



Q(x)(VS-u) = 0. (B.2) 

This leads us to define the notion of a well-behaved function on the boundary. 

Definition B.l. A function S is said to be well-behaved on the boundary provided 
Eauation \B.2\ is satisfied only at a finite number of boundary locations for all-most all 
u. 

A natural question ensues: "Why is our assumption weak?". The reason is as follows. 

Firstly note that any smooth boundary can be approximated to a given degree 
of accuracy by a finite number of hyper-planes. Hence to streamline our discussion 
here we assume that the boundary T is composed of a sequence of hyper-planes. On 
any given hyperplane Q(x) remains fixed. Recall that u ^ C, hence VS ^ u on T. 
Since rank Q is d — 1 and VS — u is required to be orthogonal to all the d — 1 rows 
of Q for condition IB. 21 to hold, VS — u is confined to an 1-D subspace. So if enforce 
VS* to vary smoothly on the hyperplane and not be stagnant, we can circumvent the 
occurrence of infinite number of second kind stationary points for all u. But as we 
mentioned above, the constraint is required only for all-most all u and not Vu. This 
facilitates us to further relax the condition. 

Let T> denote a portion of T where VS* = t for some constant vector t. Let 
u = Uo and u = Ui result in infinite number of second stationary points on T>. As 
VS — u is limited to a 1-D subspace, we must have t — Ui = A(t — uq) for some A, 
i.e, Ui = (1 — A)t + Au . So in any given region of T, there is utmost a 1-D subspace 
(measure zero) of u which results in an infinite number of second kind stationary 
points in that region. Our well-behaved condition is then equivalent to assuming that 
"the number of planar regions on the boundary where VS = constant is finite". 

To press the point home, we further illustrate it with a 2D example. Consider 
a line segment on the boundary xi = mx\ + 6. Without loss of generality assume 

1 

m 

Si + mS-2 = ui + mu2 where Si = J^. So if we plot the sum Si + 771S2 for points 
along the line, the requirement reduces to the function Si +mS2 not oscillating infinite 
number of times around infinite number of ordinate locations u\ + mui- Easy to see 
that the imposed condition is indeed weak and is satisfied by almost all functions 
barring a few exceptions. 



the parametrization y = x%. Then Q = 



Equation IB. 21 can be interpreted as 



